addpath(genpath('C:\Users\amir\Desktop\optimization_codes'))
cd C:\Users\amir\Desktop\optimization_codes
load objmat; 

[vertex,face] = read_vtk(S.tumor_type);  %read vtk image in the form of point cloud x y z
vertex=vertex';
face=face';
img1=s2v(vertex,face);
img=zeros(size(img1))+img1;
for k=1:size(img1,3)
    for j=1:size(img1,2)
        inters=find(img1(:,j,k)==1);
        [row2 col2 I32] = ind2sub(size(img1),inters);
        if(size(row2,1)>1)
            img(row2(1):row2(end),j,k)=1;
%             for i= 0:floor(size(row2)/2)-1
%                 img(row2(2*i+1):row2(2*i+2),j,k)=1;
%             end
        end
    end
end
qtumor=find(img(:,:,:)==1);
[row col I3] = ind2sub(size(img), qtumor);
exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\test_tumor1.obj', [row(:)+min(S.x) col(:)+min(S.y) I3(:)+min(S.z)], []);
if(~isempty(S.vessels))
    [vertex1,face1] = read_vtk(S.vessels);  %read vtk image in the form of point cloud x y z
    vertex1=vertex1';
    face1=face1';
    vessels=s2v(vertex1,face1);
    for k=1:size(vessels,3)
        for j=1:size(vessels,2)
            inters=find(vessels(:,j,k)==1);
            [row1 col1 I31] = ind2sub(size(vessels),inters);
            if(size(row1,1)>1)
                vessels(row1(1):row1(end),j,k)=1;
%                 for i= 0:floor(size(row1)/2)-1                 inja bayad
%                doros konam k voxelization vaghti chandta noghte mojaver
%                roshan hastan kar kone
%                     vessels(row1(2*i+1):row1(2*i+2),j,k)=1;
%                 end
            end
        end
    end
    % this part will crop artery image to 70*70*70 cube and add points if
    % outside the domain
    crop_vessellowx=round(min(vertex(:,1))-min(vertex1(:,1))-(70-(max(vertex(:,1))-min(vertex(:,1))))/2);
    crop_vessellowy=round(min(vertex(:,2))-min(vertex1(:,2))-(70-(max(vertex(:,2))-min(vertex(:,2))))/2);
    crop_vessellowz=round(min(vertex(:,3))-min(vertex1(:,3))-(70-(max(vertex(:,3))-min(vertex(:,3))))/2);
    crop_vesselupx=round(max(vertex1(:,1))-max(vertex(:,1))-(70-(max(vertex(:,1))-min(vertex(:,1))))/2);
    crop_vesselupy=round(max(vertex1(:,2))-max(vertex(:,2))-(70-(max(vertex(:,2))-min(vertex(:,2))))/2);
    crop_vesselupz=round(max(vertex1(:,3))-max(vertex(:,3))-(70-(max(vertex(:,3))-min(vertex(:,3))))/2);
    qlow=0;plow=0;tlow=0;qup=0;tup=0;pup=0;q_up=0;t_up=0;p_up=0;
    if(crop_vessellowx <=0 )
        qlow=crop_vessellowx;
    end
    if(crop_vessellowy <=0 )
        tlow=crop_vessellowy;
    end
    if(crop_vessellowz <=0 )
        plow=crop_vessellowz;
    end
    if(crop_vesselupx<=0)
        qup= crop_vesselupx;
%         q_up=-(abs(qlow)+crop_vessellowx)+size(zero_matrix,1)-70;
    end
    if(crop_vesselupy<=0)
        tup= crop_vesselupy;
%         t_up=-(abs(tlow)+crop_vessellowy)+size(zero_matrix,2)-70;
    end
    if(crop_vesselupz<=0)
        pup= crop_vesselupz;
%         p_up=-(abs(plow)+crop_vessellowz)+size(zero_matrix,3)-70;
    end
    zero_matrix=zeros(size(vessels,1)+abs(qlow)+abs(qup),size(vessels,2)+abs(tlow)+abs(tup),size(vessels,3)+abs(plow)+abs(pup));
    zero_matrix(abs(qlow)+1:end-abs(qup),abs(tlow)+1:end-abs(tup),abs(plow)+1:end-abs(pup))=vessels;
    if(crop_vesselupx<=0)
        q_up=-(abs(qlow)+crop_vessellowx)+size(zero_matrix,1)-70;
    end
    if(crop_vesselupy<=0)
        t_up=-(abs(tlow)+crop_vessellowy)+size(zero_matrix,2)-70;
    end
    if(crop_vesselupz<=0)
        p_up=-(abs(plow)+crop_vessellowz)+size(zero_matrix,3)-70;
    end
    vessel=zero_matrix(abs(qlow)+1+crop_vessellowx+q_up:abs(qlow)+crop_vessellowx+70+q_up, abs(tlow)+1+crop_vessellowy+t_up:abs(tlow)+crop_vessellowy+70+t_up,...
        abs(plow)+1+crop_vessellowz+p_up:abs(plow)+crop_vessellowz+70+p_up);
    qlow=find(vessel(:,:,:)==1);
    [row_vessel col_vessel I3_vessel] = ind2sub(size(vessel), qlow);
    S.row_vessel=row_vessel; S.col_vessel=col_vessel; S.I3_vessel=I3_vessel;
    exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\test_vessel.obj', [S.row_vessel-1+min(S.x)-S.boundx S.col_vessel-1+min(S.y)-S.boundy S.I3_vessel-1+min(S.z)-S.boundz], []);
end
x=vertex(:,1);
y=vertex(:,2);
z=vertex(:,3); 
%plot3(x,y,z,'.')
S.boundx=round((70-(max(x)-min(x)))/2);S.boundy=round((70-(max(y)-min(y)))/2);S.boundz=round((70-(max(z)-min(z)))/2);
x1=linspace(min(x)-S.boundx,max(x)+S.boundx,30);     %creat a 3D point cloud Hyper cube around tumor location
y1=linspace(min(y)-S.boundy,max(y)+S.boundy,30);
z1=linspace(min(z)-S.boundz,max(z)+S.boundz,30);
[X Y Z] = ndgrid(x1,y1,z1);
xyz = [X(:),Y(:),Z(:)];

%plot3(xyz(:,1),xyz(:,2),xyz(:,3),'.')
S.X1=xyz(:,1);S.x=x;S.face=face;      
S.X2=xyz(:,2);S.y=y;
S.X3=xyz(:,3);S.z=z;
S.result=inhull([xyz(:,1),xyz(:,2),xyz(:,3)],[x,y,z]);
S.maxX=max(x);S.minX=min(x);
S.maxY=max(y);S.minY=min(y);
S.maxZ=max(z);S.minZ=min(z);
% S.maxX=max(row)+max(x);S.minX=min(x);
% S.maxY=max(col)+max(y);S.minY=min(y);
% S.maxZ=max(I3)+max(z);S.minZ=min(z);



% exportOBJ( 'C:\Users\amir\Desktop\optimization_codes\test_tumor1.obj', [row(:)+min(S.x) col(:)+min(S.y) I3(:)+min(S.z)], []);
S.tumor_img=img;
S.row=row;S.col=col;S.I3=I3;
S.count=0;
save('objmat','S')